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Abstract 

A general procedure for constructing conservative numerical integra- 
tors for time dependent partial differential equations is presented. In 

^> particular, linearly implicit methods preserving a time discretised version 

of the invariant is developed for systems of partial differential equations 
with polynomial nonlinearities. The framework is rather general and al- 
lows for an arbitrary number of dependent and independent variables with 
derivatives of any order. It is proved formally that second order conver- 

I gence is obtained. The procedure is applied to a test case and numerical 

^— i experiments are provided. 

1 Introduction 

Schemes that conserve geometric structure have been shown to be useful when 
studying the long time behaviour of dynamical systems. Such schemes are 
sometimes called geometric or structure preserving integrators |T51 [TH] . In this 
i-H paper we shall mostly be concerned with the conservation of first integrals. 

Even if a presumption in this work is that the development of new and better 
integral preserving schemes is useful, we would still like to mention some situ- 
ations where schemes with such properties are of importance. In the literature 
one finds several examples where stability of a numerical method is proved by 
directly using its conservative property, one example is the scheme developed 
for the cubic Schrodinger equation in [11] . Another application where the exact 
• • preservation of first integrals plays an important role is in the study of orbital 

. stability of soliton solutions to certain Hamiltonian partial differential equations 

^ (PDEs) as discussed by Benjamin and coauthors [TJ[2]- 

For ordinary differential equations (ODEs) it is common to devise relatively 
general frameworks for structure preservation. This is somewhat to the contrary 
of the usual practice with partial differential equations where each equation un- 
der consideration normally requires a dedicated scheme. But there exist certain 
fairly general methodologies that can be used for developing geometric schemes 
also for PDEs. For example, through space discretisation of a Hamiltonian PDE 
one may obtain a system of Hamiltonian ODEs to which a geometric integrator 
may be applied. Another approach is to formulate the PDE in multi-symplectic 
form, and then apply a scheme which preserves a discrete version of this form, 
see |3] for a review of this approach. 
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In this paper we consider methods for PDEs that are based on the discrete 
gradient method for ODEs. The discrete gradient method was perhaps first 
treated in a systematic way by Gonzalez [16] , see also [HI [25] . For PDEs one 
may derive discrete gradients either for the abstract Cauchy problem, where 
the solution at any time is considered as an element of some infinite dimen- 
sional space, or one may semidiscretise the equations in space and then derive 
the corresponding discrete gradient for the resulting ODE system. This last 
procedure has been elegantly presented in several articles by Furihata, Matsuo 
and collaborators, see e.g. [T2j [El O ED E2 E31 [24] , using the concept of dis- 
crete variational derivatives. See also the monograph [15]. The first part of this 
paper develops a similar framework that is rather general and allows for an ar- 
bitrary number of dependent and independent variables with derivatives of any 
order. The suggested approach does not require the equations to be discretised 
in space. 

We consider a class of conservative schemes which are linearly implicit. By 
linearly implicit we mean schemes which require the solution of precisely one 
linear system of equations in each time step. This is opposed to fully implicit 
schemes for which one typically applies an iterative solver that may require a 
linear system to be solved in every iteration. For standard fully implicit schemes 
one would typically balance the iteration error in solving the nonlinear system 
with the local truncation error. However, for conservative schemes the situation 
is different since exact conservation of the invariant requires that the nonlinear 
system is solved to machine precision. This work can be seen as a generalisation 
of ideas introduced in Section 6] . 

It may not, in general, be an easy task to quantify exactly what can be 
expected of gain in computational cost, if any, when replacing a fully implicit 
scheme with a linearly implicit one. For illustration we present an example 
where the KdV equation 

u t + u xxx + (u 2 ) x = (1.1) 

is solved on a periodic domain using a fully implicit scheme 
jjn+l _ jjn jjn+l + jjn^ / gjn+iy + jjn+ljjn + gjn}2 



At 2 
and a linearly implicit scheme 



= 0, (1.2) 



jjn+2 _ jjn jjn+2 + jjn / jjn+2 + jjn+l + Jjn\ 

_ + 2 + J =0, (1.3) 

where U n {x) w u(x, t n ) = u(x, t° + nAt). 

These schemes are derived in Section [6] For space discretisation centered 
differences are used for both schemes. Note that the linearly implicit scheme 



(1.3) has a multistep nature, but should not be confused with standard linear 
multistcp methods. Furihata et al. sometimes use the term "multiple points 
linearly implicit schemes" to emphasise this fact. 

The schemes are both second order in time, but in our example the linearly 
implicit multistep scheme has an error constant which is about 3-4 times larger 
than the fully implicit one-step scheme. In Figure |1.1| we plot the global error 
versus the number of linear solves for the two schemes. The linearly implicit 
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scheme solves one linear system in each time step. The fully implicit scheme, 
on the other hand, solves a linear system for each Newton iteration which is 
repeated to machine precision in each time step. For the largest time step in 
our experiment this amounts to 561 linear solves per time step. The linear 
systems in each of the two cases have the same matrix structure, they are both 
penta-diagonal, and we therefore assume that the cost of solving the linear 
system is approximately the same for both methods. The s-axis in Figure |1.1| 
can thus be interpreted as a measure of the computational cost in each scheme. 
The plot shows that for a given global error the linearly implicit scheme is 
computationally cheaper than the fully implicit scheme. 

There are situations in which the results from this example may be less 
relevant. For instance, the iteration method used in a fully implicit schemes 
may use approximate versions of the Jacobian for which faster solvers can be 
applied, therefore the cost of a linear solve may not be the same for the two 
types of schemes. For large time steps, both types of schemes are likely to en- 
counter difficulties, but for slightly different reasons. The fully implicit scheme 
may experience slow or no convergence at all of the iteration scheme, whereas 
the linearly implicit scheme may become unstable for time steps over a certain 
threshold [Hj ■ For instance, in the case of the stiff ODE considered by Gonzales 
and Simo [17] one will observe that the stability properties of the linearly im- 
plicit schemes will be completely lost whereas fully implicit conservative schemes 
behave remarkably well. For large-scale problems, one may have situations in 
which iterative linear solvers are required and where one cannot afford to solve 
these systems to machine accuracy, in such cases the linearly implicit schemes 
are less useful. In conlusion, we believe that which of the two types of schemes 
that is preferable depends on the PDE and the circumstances under which it is 
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to be solved. 

The two schemes used in the example above have slightly different conser- 



vation properties. The first one (1.2 1 conserves the exact Hamiltonian 



7/7 "]= / [!([£)* -I([7»)3) ,Lr. 



£2 



whereas the second scheme (1.3) conserves what we will define as the polarised 
Hamiltonian 



| ((f/") 2 + {u™ +1 ) 2 ) ~ i ({U n ) 2 U n+1 + (U n+1 ) 2 U n )^j 



H[u n , u n+L ] = j ^- {(uzy + (u% +i y) - - {{u n yu n+L + (u n+L yu n ) ) &x. 

Both of these functions are approximations to the true Hamiltonian, the first is a 
spatial approximation for a fixed time, and the second also includes an averaging 
over time. The intention is that in both cases one can see the methods as exactly 
preserving a slightly perturbed first integral over very long times. That this 
seems to work for the chosen example is clearly seen in the first plot in Figure 



1.2 where we plot the error in H as a function of the solution obtained by the 



linearly implicit scheme (1.3 1. We integrate to t = 1000 and the error is plotted 
from t = 980 to t — 1000. Notice that in this example there is no drift in the 
energy error. The corresponding error plots for T-L as a function of the solution 



of (1.2) and H as a function of the solution of (1.3) are omitted since they are 



both preserved up to round-off error, and thus not that interesting. The second 



plot in Figure 1.2 shows how the error in H at the endpoint depends on At. 



Empirically, we have the relation 

U[U n ] =H[U°] + C(At) 2 , 

where C is a constant that depends on the solution, but not on n. See Section 
[6] for another example that tests the long time structure preserving properties 
of these schemes. 

This and similar examples show that there are situations where linearly 
implicit schemes can be a better choice than their fully implicit counterparts. 
Figure |1.1| shows that the linearly implicit scheme is cheaper, while Figure |1.2| 
shows that both solutions have similar long-term behaviour. Similar favourable 
behaviour of linearly implicit schemes can be found in the literature. For the 
cubic Schrodinger equation there are such conservative schemes based on time 
averaged versions of the Hamiltonian by Fei et al. [11] and by Besse [3]. Exam- 
ples of methods for other PDEs can be found in the monograph [TS] and the 
papers gU] and gS]- 

In the next section we define the PDE framework that we use. Then, in 
Section [3] we consider discrete gradient methods and how they can be applied 
to PDEs. We study in particular the average vector field method by Quispel 
and McLaren [STJ and the discrete variational derivative method by Furihata, 
Matsuo and coauthors [H EH HQ 121 121 121] ■ We develop a framework that 
works for a rather general class of equations. 

The key tools for developing linearly implicit methods for polynomial Hamil- 
tonians are treated in Section [1J introducing the concept of polarisation. There 
is some freedom in this procedure, and we show through a rather general exam- 
ple term how the choice may significantly affect the stability of the scheme. 
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Figure 1.2: The error in % as a function of the solution obtained by the linearly 
implicit scheme 1.3 as a function of t (top) and At (bottom). The dotted line 
is a reference line C(At) 2 . 



We defer the introduction of spatial discretisation until Section [5] This is 
done mostly in order to keep a simpler notation, but also because our approach 
concerns conservative time discretisations and is essentially independent of the 
choice of spatial discretisation. The last section offers some more details on 
the procedure for constructing schemes and we give some indication through 
numerical tests on the long term behaviour of the schemes. 



2 Notation and preliminaries 

We consider integral preserving PDEs written in the form 

X'TJ 

u t =V~, (2.1) 
ou 

n[u]= [ g[u]dx= [ g((u a j))dx, ncu d , (2.2) 

is the preserved quantity and T> is a skew-symmetric operator that may depend 
on u. We write dx = dx 1 ■ ■ ■ dx d . We remark in passing that the class of PDEs 



where 



which can be written in the form (2.1 1 contains the class of Hamiltonian PDEs, 
however we do not make the additional assumption that T> satisfies the Jacobi 
identity [55]. By (it") we mean u itself, which may be a vector u — (u a ) £ 
M. m , and all its partial derivatives with respect to all independent variables, 
[x 1 , . . . ,x d ), up to and including some degree v. Thus, J is a multi-index, we 
let J = (Ji, . . . , j r ), where r = \ J\ the number of components in J, and 

d r u a 



dx^ 1 ■ ■ • dxi r 



< r < v. 
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As in [25], the square brackets in (2.2 1 are used to indicate that a function 



depends also on the derivatives of its arguments with respect to the independent 
variables. In one dimension d — 1 and m — 1, for example, we have 

ow-o ««,»-<>(., £,...,§£ 

The variational derivative is an m- vector depending on u" for | J| < v' where 
v' > v. It may be defined through the relation [25J p. 245] 



SH , a 

— • <p Ax = — 
a ou de 



H[u + etp], (2.3) 



c=0 



for any sufficiently smooth m- vector of functions <p(x). One may calculate 
by applying the Euler operator to Q[u], the a-component is given as 



5H 

8u 



E a G[ U ], (2.4) 



where 

E Q = £ (-l) |J '^i (2-5) 

OUj 

\J\<U J 

so that the sum ranges over all J corresponding to derivatives u" featuring in 
Q. We have used total derivative operators, 

it, n , ...n,. a = V?^- 

a, J J 

In parts of the paper we refer to Hamiltonians as polynomial, or specifically 
quadratic. By this we mean that % is of a form such that Q is a multivariate 
polynomial in the indeterminates uj, which in the quadratic case is of degree at 



most two. For example, the KdV equation (1.1 1 has a polynomial Hamiltonian 
of degree 3 

K\n\ : [ ) dr. 



In this case Q = 0(u, u x ) and thus m = d = v = 1, and we get 

m writ ^ dG d dG a*\ 
— =Eg((u >)) = ---— (2.6) 

= -u 2 - u xx . (2.7) 

We always assume sufficient regularity in the solution and that the boundary 
conditions on D, are such that the boundary terms vanish when doing integration 
by parts, for example periodic boundary conditions. The operator T> should be 
skew-symmetric with respect to the L 2 inner product 



(Vv)wdx = - / v(Vw)dx Vu,w. (2.8) 
For the KdV case we simply have T> = 
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Furthermore, to be a true Hamiltonian system it should induce a Poisson 
bracket on the space of functionals as described e.g. in [23 Ch 7.1], meaning 
that the Jacobi identity must be satisfied. However, the approach presented here 
only requires T> to be skew-symmetric so that the functional % is a conserved 
quantity. In the case that the PDE has more than one Hamiltonian formulation, 
we may make a choice of which of the integrals to preserve. Our approach does 
not in general allow for the preservation of more than one Hamiltonian at the 
same time, for this see the upcoming paper [lOj . 

PDEs such as the wave equation are typically written with u u appearing on 
the left hand side, in such cases we double the dimension of u in order to apply 
the stated framework. For complex equations one may do something similar, 
splitting either into a real and an imaginary part, or adding in the complex 
conjugate as a separate variable. 

3 Discrete gradient and variational derivative 
methods 

Discrete gradient methods for ODEs were introduced by Gonzalez [Tfo] . See also 
[6], [I], [25], and [HI Chapter V.5]. Recently this idea has been applied to PDEs 
in the form of the average vector field (AVF) method [3] and in a somewhat 
more general setting, the discrete variational derivative (DVD) method. 

We recall the definition of a discrete gradient as presented for ODEs. If 
H : R M -» R, a discrete gradient is a continuous map V : R M x R M R M 
such that for every u and v in R M 

H(u) - H(v) = Vi? (v, u) • (u - v), 
Vff(u,u) = V/f(u). 

Since an ODE system preserving H can be written in the form 

^ = S(y)Vff(y) 

for some skew-symmetric matrix S(y), one obtains a conservative method simply 
by defining approximations y™ ~ y(t n ) = y(t + nAt) through the formula 

n+l _ n 

At = 5Vg (y"'y )» 

where 5 1 , typically allowed to depend on y™ and y n , is some skew-symmetric 
matrix approximating the original S. 

There are many possible choices of discrete gradients for a function H , see 
for instance [T51 [5S] . A particular example is the one used in the AVF method 
defined as ^ 

V AVF ff(v,u) = f V#(£u+(l-£)v)de 
Jo 

When applying this approach to PDEs the obvious strategy is to discretise the 
Hamiltonian H[u\ in space, replacing each derivative by a suitable approximation 
like e.g. finite differences, to obtain a Hamiltonian "Hd(u) as for ODEs. Similarly, 
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the skew-symmetric operator V is replaced by a skew-symmetric M x M-matrix 
T>d to yield the scheme 



— = V d Vft d (u", u" +i ) (3.1) 

for advancing the numerical solution u n at time f™ to u n+1 at time t n+1 . Ex- 
amples are worked out for several PDEs in [3] . 

Furihata, Matsuo and coauthors present a whole framework for discretis- 
ing PDEs in the variational setting in a series of papers, providing a discrete 
analogue of the continuous calculus, see for instance [T2]. They discretise Q to 
obtain Qd using difference operators, and then the integral in % is approximated 
by a sum to yield "H^. Then they derive a discrete counterpart to the variational 
derivative, and finally state the difference scheme in a form which is a perfect 
analogue to the Hamiltonian PDE system (2.1), letting 

u» +1 - u" 6H d 



At ^(u",^ 1 )' 

The use of integration by parts in deriving the Euler operator is mimicked 
by similar summation by part formulas for the discrete case. Their discrete 
variational derivative is in fact rather similar to a discrete gradient, as it satisfies 
the relation 

H d (u)--H d (v) = {J^,u-v) (3.2) 
d(v,u) 

for the discrete L 2 inner product. 

In the present paper, we focus on the time dimension in most of what follows, 
thus we shall defer the steps in which % and thereby Q are discretised in space. 
But (3.2) makes perfect sense after removing the subscript d, replacing u and v 
by functions u and v, and the discrete L 2 inner product by the continuous one. 
A discrete variational derivative (DVD) is here defined to be any continuous 
function j^^y of (u!"',f') satisfying 

f 81-1 

H[u] - U[v) = /- Au-v)dx, (3.3) 

Jq(>{v,u) 

(3.4) 



5H 5H 



5{u, u) Su 

The integrator yields a continuous function U n := U n (x) « u(x,t n ) for each t n 

U n+1 - U n _ SH . . 

At ~ <5(/7", [/"+ 1 )- [ j 



By combining (3.3) and (3.5) we see that the method preserves T~L. 

The AVF scheme can of course also be interpreted as a discrete variational 
derivative method where 

,m -< r ™[tu+(l-Qv]dt. (3.6) 



S(v,u) J Q Su 



The fact that (3.6) verifies the condition (3.3) is seen from the elementary iden- 
tity ^ 

H[u] -H[v] = J ±H[e u + (i-£) v ]6Z. (3.7) 
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The derivative under the integral is written 



H[v + + e)(u - v)] 

6=0 

— [f u + (1 - f )«] (« - w) dx. 
ou 



Now substitute this into (3.7) and interchange the integrals to obtain (3.3 1. 

In most of the cited papers by Furihata, Matsuo and coauthors, the notion 
of a DVD method is less general than what we just presented, in the sense that 



the relation (3.2 1 is not actually used as the defining equation for a discrete 
variational derivative. Instead the authors present a relatively general format 
that can be used for discretising H, this format is depending on the class of PDEs 
under consideration, and they work out the explicit expression for a particular 
discrete variational derivative. To give an idea of how the format may look 



like, we briefly review some points from [TJ] where PDEs of the form (2.2) are 
considered with d = m = v = 1 such that Q = Q(u,u x ). Q is assumed to be 
written as a finite sum 

G(u, u x ) = ^2 aefe(u)ge(u x ). (3.8) 



where fg and g% are differentiable functions. Q The form (3.3) is then derived 
through 

t< \ , ^ ft\(\ M u ) - f*( v ) 9e(u x ) + 9t{v x ) , , 

ft{u)gi{u x ) - fe{v)g e \v x ) = (it - v) 

u — v 2 

. 9e(u x ) - gt{v x )ft{u) + fe(v) . , 
H -z \ u x - V x) 

followed by an integration by part on the second term. This technique can be 
extended in any number of ways to allow for more general classes of PDEs. For 
instance, one may allow for more factors in (|3.8|), like 



G[u] = ^2 a t Y[ 9e,j{dju) 



and repeated application of the formula ab~cd = (b — d) + (a — c) to this 
equation combined with integration by parts will result in a discrete variational 
derivative. 

Schemes which are built on this particular type of discrete variational deriva- 
tive will be called the Furihata methods in the sequel since it was first intro- 
duced in [12] . Matsuo et al. extend the method to complex equations in [22], 
while [131 121) derive methods for equations with second order time derivatives. 
Other papers using the discrete variational derivative approach include [14] . 
[23, and 0D]. 

The lack of a general formalism in the papers just mentioned, makes it 
somewhat difficult to compare the approach to the AVF method and characterise 



In 1121 the expression is discretised in space and gi(u x ) is replaced by a product 
9l (^fc~^fc)9<r (^k ^ fc ) wnere Su a- n d <5j7 are forward and backward divided differences respec- 
tively. 
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in which cases they lead to the same scheme. Taking for instance the KdV 
equation one easily finds that both approaches lead to the scheme (1.2 1, 
however, considering for instance the Hamiltonian 

% [u] = / tiitj dx 

JO 

one would obtain two different types of discrete variational derivative in the Fu- 
rihata method and the AVF method, that is 7^ TCfaf • ^ n some important 

cases, the Furihata method and the AVF method lead to the same scheme. 

Theorem 3.1. Suppose that the Hamiltonian H[u] is a linear combination of 
terms of either of the types 

1. J n dju ■ &Kudx for multi-indices J and K, or 

2. Jq g{dju) dx for differ entiable g : K — > K. 

Then the AVF and the Furihata methods yield the same scheme 
Proof. It suffices to check one general term of each type. 



1. We find the variational derivative using (|2.3|) 

sn 

Su 



((-l)l'l + {-l)W)d J+K u. 



Inserting the variational derivative into (3.6) gives 

5H AVF ( , _ _ 1 _ri . / - % 1 k\ \ ^ f u + v 



8(v, u) 



((-!)'•" + (-1)'*') 9,+K 



To find the discrete variational derivative of the Furihata method we com- 
pute 

W[u] — J H[v]= / dju ■ &ku — djv ■ dxv dx 
Jn 

z Jn 

After integration by parts we get 

o(v,u) V / 

and we see that — r = — . 

o(v,u) o{v,u) 

2. In this case we get 



so that 

5U f 



S(v, u) 



f = (-ir\d jg '(d JU ), 

ou 



(-1)1 J l / djg'(dj(Zu+(l-Ov))dt 



g(d.;u) - g(djv) 
dju - djv 
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□ 



For the Furihata method one would here just compute 

H[u] ~ H[v\ = / — (dju - djv) dx 

J n dju - djv 

and integration by parts yields ^ F = ^ AVF , ■ 

o(v,u) o(v,u) 

4 Linearly Implicit Difference Schemes 
4.1 Polarisation 

The key to constructing conservative linearly implicit schemes will be to portion 
out the nonlincarity over consecutive time steps. In effect, this means that we 
replace the original Hamiltonian % with an approximate one H . We shall call 
H a polarisation of % since its definition resembles the way an inner product is 
derived from a quadratic form. We shall see that the difference scheme resulting 
from such a polarised Hamiltonian will be a multistep method. This method 
will now preserve exactly H, as opposed to % for the methods in the previous 
section. The requirements on H are given in the following definition. 

Definition 4.1 (The polarised Hamiltonian). Given a Hamiltonian T-L[u] the 
polarised Hamiltonian H depends on k arguments, and is: 

• Consistent 

H[u,u,...,u] =n[u}. (4.1) 

• Invariant under any cyclic permutation of the arguments 

H[wi,w 2l ...,w k ) = H[w 2 , . . .,w k ,Wi]. (4.2) 
Polarisations exist for any Hamiltonian, this is asserted by the example 
H[w 1> w 2 ,...,w k ] = -{n[w 1 ]+n[w 2 ] + ---+H[w k }). 
We may impose the polarisation directly on the density Q((v,j)), letting 

H[wi,w 2 ,...,w k ]= / G[w 1 ,w 2 , ■ ■ ■ ,w k ]dx. 



The conditions (4.1), (4.2) are then inherited as 



G(u, u,...,u) = G(u), G[w 1 ,w 2 , ...,w k ] = G[w 2 , . . . , w k> wi] 



In Section [4~2] we will discuss local order of consistency, it will then be conve- 
nient to make the stronger assumption that % and H are at least twice Frechet 
diffcrentiable. To distinguish from the weaker notion of variational (Gateaux) 
derivative, we replace S by d, noting that the first derivative in the two defi- 



nitions are the same when they both exist. We then find from (4.1) and (4.2) 
that the Frechet derivatives satisfy the relation 

^t»>=W" * (43) 
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For the second derivatives, we find the identity 
d 2 H , , d 2 H 



[U,. 



dwidwj ' ' dwidwk+2-j 
which is used to compute 



. , [k/2\ + 1, (4.4) 



dv? 



u = < 



d 2 H 


k+l 
2 

1=2 


d 2 H 


dw\ 


dwidwi 


d 2 H 


k 
2 

e=2 


d 2 H 


dwf 


dwidwi 



+ 



d 2 H 



k odd, 



k even, 



(4.5) 



all second derivatives on the right being evaluated at [it, . . . , u\. 
4.1.1 Polynomial Hamiltonians 

The polarisation of polynomial Hamiltonians will be key to constructing linearly 
implicit schemes. We will now explain in detail how to do this, and we begin 
with an example term in the integrand Q[u\ = G{dju) depending on just one 
scalar indeterminate, namely G{z) = z p where z = dju a for some (J, a) and 
where p < 4. This example is important not only as a simple illustration of 
the procedure, but also because terms of this type are common in many of the 
Hamiltonians found in physics. As we will see in the next section, it will be 
natural to use two arguments, k — 2, in the polarised Hamiltonian. In fact, we 
need to restrict ourself to cases with polynomial Hamiltonians for our technique 
to yield linearly implicit schemes. Then, by using k > [p/2], we can obtain 
polarised Hamiltonians which are at most quadratic in each argument. We call 
these quadratic polarisations. We see that if k = 2 then cyclic is the same 
as symmetric G(u,v) — G{v,u) 1 and the possible quadratic polarisations for 
p = 2, 3, 4 are respectively, 



p = 2 

p = 3 
p — 4 



G(u, v) = 



, u 2 + v 2 



G(u, v) = uv 



2 

U + V 



+ (l-6)uv, 9 e [0,1], 



G(u, v) — u 2 v 2 



(4.6) 

(4.7) 
(4.8) 



Note that for these monomials both the third and fourth degree case are uniquely 
given, but the second degree case is not. In Section |4.3| we will consider how 
the choice of 9 influences the stability of the scheme for a term which appears 
frequently in PDEs. 

We now consider the general case when Q[u] is a multivariate polynomial 
in N v variables of degree p. It suffices in fact to let G[{uj)) be a monomial 
since each term can be treated separately, for u £ M. Nv . For a convenient 
notation, we rename the vector of indeterminates (uj) by using a single index 
i.e. u = (iti, . . . , un u ) and write 

G{u) = u Zl u i2 . . .u ip . (4.9) 

One may use the following procedure for obtaining a quadratic polarisation from 
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1. Group the factors of the right hand side of (4.9| into pairs z r = Ui 2r _ 1 Ui 2r 
and if p is odd z k — Ui p . Set 

K(zi, ...,z k ) = z x ---z k . 



Note that there are potentially many ways of ordering the factors in (4.9) 
which give rise to different polarisations. 

2. Symmetrise K with respect to the cyclic subgroup of permutations. Let- 
ting the left shift permutation a be defined through crK{z\, . . . , z k ) — 
K(z2,...,Zk,zi), we set 

1 k 

G(z u ...,z k ) = rYj ak 1 K(z 1 ,...,z k ). 

fe=i 



The resulting G is now both consistent (4.1) and cyclic (4.2) 



4.2 Linearly Implicit Methods 

We may now define the discrete variational derivative for this polarised Hamil- 



tonian as a generalisation of (3.3) and (3.4). We let 



SH 



S(w x , . . . , w k+ i) 
be a continuous function of k + 1 arguments, satisfying 

f SH 

H[w2,--.,W k+1 }-H[w 1 ,...,W k )= / — r 

J u d{w 1 ,...,w k+1 ) 
SH 8H 



(w k+ i -w^dx, (4.10) 



(4.11) 



S(u, . . . , u) 5u 

Our standard example will be a generalisation of the AVF discrete variational 
derivative, which we define as 



^AVF 

6(wi, . . . , w k+1 ) 



1 SH 
) Jo Swi 



[£,w k+1 + (1 - g)w l ,W2 ) ...,w k ]d£. 



(4.12) 



Here the variational derivative on the right hand side, 4^- is defined as before, 
considering H as a function of its first argument only, leaving the others fixed. 
Similar discrete variational derivatives could be derived in a number of different 
ways. In particular one finds that when the function H is quadratic in all its 
arguments, the approach used in deriving the Furihata methods would lead to 
a discrete variational derivative which is identical to that of the AVF-method. 

Now we define the polarised discrete variational derivative scheme and prove 
that, under some assumptions, this scheme is conservative, linearly implicit and 
has formal order of consistency two. 



Definition 4.2. For a Hamiltonian PDE of the form (2.1 1, let H be a polarised 



Hamiltonian of k arguments, satisfying (4.1) and (4.2), and suppose that ap- 
proximations U 3 to u(jAt, •) are given for j = 0, . . . , k — 1. 
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The polarised DVD (PDVD) scheme is given as 



U 



n-\-k 



u n 



kAt 



= kD 



SH 



5{u n , . . . ,u n + k y 



n > 0. 



(4.13) 



• D is a skew- symmetric operator approximating T>. In (4.13), D may de- 
pend on [7 n +J , l<j<fc — 10 and be consistent 

D[u,...,u] =V[u\. (4.14) 

D is called cyclic if 

D[wi,w 2 , ■ • • ,w k -i] = D[w 2 t ■ ■ ,w k -i,wi\. (4-15) 



• If the discrete variational derivative is given by (4.12 1, then the scheme is 
called the polarised AVF (PAVF) scheme. 



Theorem 4.3. The scheme (4.13) is conservative in the sense that 

rfcl 



H[U n+ \...,U n+k ] = H[U ,.. 
for any polarised Hamiltonian function H . 



Vn > 1. 



Proof. By induction, this is an immediate consequence of (4.10) 



□ 



In a framework as general as the one presented here, it is not possible to 
present a general analysis for convergence or the order of the truncation error. 
However, it seems plausible that a necessary condition to obtain a prescribed 
order of convergence can be derived through a formal Taylor expansion of the 
local truncation error, we denote this the formal order of consistency. 

Theorem 4.4. 

• The PAVF scheme has formal order of consistency one for any polarised 
Hamiltonian, and skew- symmetric operator D satisfying (4.14). 



If in addition (4.15) is satisfied, the scheme has formal order of consistency 
two. 



Proof. We show that when the exact solution is substituted into (4.13) where 
the discrete variational derivative is given by (4.12 ), then the residual is 0(At 2 ). 
Throughout the proof we assume the existence of Frechet derivatives. Writing, 
for any j, vP = u(-,P) for the exact local solution at t = P , we get for the left 
hand side 



kAt 



0{At 2 ) = V 



&H 
du 



/.■A/ fdv^^m 

du du 



+ 0(At 2 ). (4.16) 



2 D should not depend on (/ n + fe since otherwise the method would no longer be linearly 
implicit. 



14 



Next we expand (4.121 to get 



SH A 



,n+k\ 



5{u" 



where u = (u n , . . 

SH AVF 



1 SH 

Swi 



OH 




i 


d 2 H 




+ 

U 







(£u n+k + (l-£)u n ,u n+1 
d 2 H 



,n+k-l 



3=2 



dwidwj 




,u n ). Using first (4.4) and then (4.3), (4.5) we find 



S(u n , . . 



,n-\-k 



i m 

k du 



At cPH 
T Ihi 2 



(;d t U n ) + 0(At 2 ). 



(4.17) 



Expanding D we get 

u n+*-l] = p[ u n] 



D\u 



n+l 



At 



5^4 dD 

3 = 1 



8wa 



(d t u n ) + 0(At 2 ). (4.18) 



If the cyclicity condition (4.15) holds for D, we can simplify (4.18) to obtain 



D\u 



rt+l 



.,U 



n+fe-ll _ 



V\u 



k{k-l)At dD 

~2~ ~du 



(d t u n ) + 0(At 2 ) 
{d t u n )+0(At 2 ). (4.19) 



By substituting into (4.13) the expressions (4.16), (4.17) and (4.19), all terms 
of zeroth and first order cancel and we are left with 0{(At 2 )). 



□ 



Theorem 4.5. Suppose that the polarised Hamiltonian H is a quadratic poly- 
nomial in each of its arguments, then the PAVF scheme is linearly implicit 



Proof. Since H is at most quadratic in the first argument, it follows from (2.5) 

SH 
8wi 

SH 



that is of degree at most 1 in its first argument, and so we see from (|4.12|) 
that 



is linear in U 



n+k 



S(U n ,...,U n + k ) 
Since D does not depend on JJ n+k we conclude that the 



scheme (4.13) is linearly implicit. 



□ 



In some cases one wishes to have time-symmetric numerical schemes, see for 
example |18j . The numerical scheme (4.13) will in general not be symmetric, 
however it is not hard to modify the procedure to yield symmetric schemes. 
One needs to polarise H such that H is invariant also when the order of its 
arguments is reversed, it turns out that this can be achieved by symmetrising 
over the dihedral group rather than just the cyclic one. A similar adjustment 
must be made for D. 

We remark that one can construct explicit schemes by using p time steps (as 
opposed to k) in H such that H becomes p-linear (as opposed to fc-quadratic) . 
The rest of the procedure for the explicit case is the same as for the linearly 
implicit case. Clearly, one expects that explicit schemes will have more severe 
stability restrictions than the linearly implicit ones. 
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Since these multistep schemes need the k previous values, it is not self- 
starting. We have to provide the starting- values U 1 , . . . , U k in addition to the 
initial value U°. Usually these are computed using another sufficiently accurate 
conservative scheme, such as for example the AVF scheme. Another possibility 
is to use any integrator and integrate to machine precision. 



4.3 Stability 

In [5] we studied linearly implicit schemes for the cubic Schrodinger equation, 
and found that two-step schemes can develop a two-periodic instability in time. 
We also saw that this can be remedied by choosing a different polarisation of 
the Hamiltonian. 

As it turns out, a common case is when the Hamiltonian is a univariate poly- 
nomial of degree 4 or less. If we polarise this Hamiltonian using two time-steps, 



we get three linearly independent H, corresponding to (4.6), (4.7), and (4 
The third and fourth degree Hamiltonians are uniquely given. However, in the 
second degree case we can choose 9 € [0, 1] such that the scheme becomes un- 
stable. Since Hamiltonians of the type ( |4.6[ ) appear in many important PDEs 
it may be useful to determine which 9 6 [0, 1] lead to unstable schemes. 
We choose to study the test equation with Hamiltonian 

H[u] = ^ J^u 2 x dx, 

and a skew-symmetric operator T> which satisfies the eigenvalue equation 

Ve lkx = iX k e ikx , X k e R 
for all integers k. The Airy equation 

u t + u xxx = 

is of this type with T> = — d x and X k = —k. Other equations which have such 
terms in the Hamiltonian include the nonlinear Schrodinger equation, the linear 
wave equation, the KdV equation, and the Kadomtsev-Petviashvili equation. 



Rewriting (4.6) gives 



H[v, u] = \j + (1 - e K«*) 

And the numerical scheme is 

Tjn+2 _ Tjn f jjn+2 , jjn \ 

2At = -V [e - 2 + ** + (1 - ffjU^j . (4.20) 

Since this is a linear equation we can use von Neumann stability analysis [S]. 
We insert the ansatz 

U n {x) = Ce ikx 
to obtain the quadratic equation 

(1 - 9t\)C, 2 - 2(1 - e)ri(- (1 + 9tV) = 0, r = \ k Atk 2 . (4.21) 
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-6 -4 -2 2 4 6 



Figure 4.1: The numerical solution of the Airy equation with two different values 
of 9. The two solutions are shown after n = 10 6 time steps (9 = 0.5) and n = 115 
time steps {9 = 0.49). 



A necessary condition for stability is \(\ < 1 which implies 



1 1 

" 2 ~ 2^2' 



Assuming that {\kk 2 }kez is unbounded, we must require that 9 is chosen 
greater than or equal \. This is exactly the condition found in [9] for the cubic 



Schrodinger equation. When 9 > \ the roots of (4.21) satisfy = \Qi\ = 1- 
In Figure 4.1 we solve the Airy equation with the scheme ( 4.20[ ) using 9 = 0.5 



and 9 = 0.49. We use the initial value u(x, 0) = sin(x), which, in the exact case, 
yields the traveling wave solution u(x,t) — sin(x + t). The 9 = 0.49 solution 
blows up in few time steps, while the 9 — 0.50 solution shows no signs of 
instability. Doing a discrete Fourier transform of the unstable solution we see 
that the instability starts at high frequencies, that is large k, which corresponds 
to the results shown above. 

There might be cases where the scheme develop instabilities due to spurious 
modes no matter what polarisation one chooses. A full stability analysis of 
either the fully or linearly implicit schemes is to our knowledge not been done. 
Standard linearisation techniques will usually lead to the conclusion that the 
schemes are neutrally stable. The nonlinear effects, however small, may still 
cause the scheme to be unstable. The tests we have done on a wide range of 
PDEs seem to indicate that the stability usually is very good, future work may 
shed a light on this issue. 
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5 Space discretisation 



Until now we have mostly considered the situation where the PDE is discretised 
in time while remaining continuous in space. The methodology developed in 
the previous sections apply equally well to systems of ODEs. Arguably, the 
most straightforward approach is simply to discretise the space derivatives in 
the Hamiltonian, for instance by finite differences. This leads to 

H(u)—>%i(u). 

One also needs to replace the skew-symmetric operator I? by a skew-symmetric 
matrix T>d- The fully implicit method reviewed in Section [3] is then just the 



discrete gradient method (3.1), which conserves the discretised Hamiltonian 
Hd{u) in every time step. 

We consider now finite difference approximations. The function space to 
which the solution u belongs, is replaced by a finite dimensional space with 
functions on a grid indexed by I g C Z d . We use boldface symbols for these 
functions. Let there be N r grid points in the space direction r so that N = 
Ni ■ ■ ■ N d is the total number of grid points. We denote by u Q the approximation 
to u a on such a grid, and by u the vector consisting of (u 1 , . . . , u m ). We will 
replace each derivative u" by a finite difference approximation Sju a , and replace 
the integral by a quadrature rule. We then let 

H d (u) = J2UQd((Sju))) i Ax. (5.1) 
iei g 

Here Aa; is the volume (length, area) of a grid cell and b = are the 

weights in the quadrature rule. The discretised Qd has the same number of 
arguments as Q, and each input argument as well as the output are vectors in 
K. . We have here approximated the function u j by a difference approximation 
Sju a , where 5j : R N — > R N is a linear map. As in the continuous case, we 
use square brackets, say F[u], as shorthand for a list of arguments involving 
difference operators F[u) — F(u, S^u, . . . , 6j u). We compute 

H d [u]-n d [v} = 

E bi E Si (^) , K u + " 0v]d£(fc(u" - v«)) Ax 



where 



, 5H d 
M(v,u) 

fiy=E^<f||^ + (i-eK]d, 



,u-v) (5.2) 



J,a 



B is the diagonal linear map B — diag(&i), i £ I g , and the discrete inner product 



used in (5.2 1 is 

(u,v)= £ ufvf. 



a.iel„ 



Notice the resemblance between the operator acting on Qd in (5.2) and the 



continuous Euler operator in (2.5). Alternatively, suppose that 



18 



1. The spatially continuous method (3.5| (using (3.6)) is discretised in space, 
using a skew-symmetric T>d and a selected set of difference quotients Sj 
for each derivative dj. 



2. Considering (2.4 1 and (2.5), the choice of discretisation operators Sj used 
in dG/duj[u] is arbitrary, but the corresponding Dj is replaced by the 
transpose <5j. 

In this case, using the same T>d, an identical set of difference operators in dis- 
cretising % (5.1), and choosing all the quadrature weights b[ = 1 the resulting 
scheme would be the same as that given by procedure outlined in the two points 
above. That is, one can get the same scheme by either discretising the Hamil- 
tonian in space first (and then deriving the scheme) or discretising the scheme 
in space first (and then deriving the conserved Hamiltonian) . 

Letting the rth canonical unit vector in M. d be denoted e r , we define the 
most used first order difference operators 



(a+u), = 

(<5«u)i = 



u 



i+e r 



Uj 



Ax r 

Uj Ui 



Ax r 



2Aa 



These difference operators are all commuting, but only the last one is skew- 
symmetric. However, for the first two one has the useful identities 



(*+) T = s; 



(s-f = -s: 



Higher order difference operators Sj can generally be defined by taking compo- 
sitions of these operators, in particular we shall consider examples in the next 
section using the second and third derivative approximations 

We may now introduce numerical approximations U™ representing the fully 
discretised system, the scheme is 



U 



Tl + l 



U" 



At 



V, 



5H d 



5(U",U"+ 1 )' 



The conservative schemes based on polarisation are adapted in a straight- 
forward manner, introducing a function Hd[wi, ■ ■ ■ , wj which is consistent and 
cyclic as in (4.1), (4.2), and a skew-symmetric map Dd depending on at most 
k — 1 arguments. The scheme is then 



Tjn+fc _ Tjr, 



kDd 



kAt '"~ a 6(U n ,...,U n + k ) 
This scheme conserves the function H d in the sense that 
H d [U n+ \ . . . ,U" +fc ] =H d [U ,..., U fc_1 ], 



(5.3) 



n > 0. 
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6 Examples 



To illustrate the procedures for constructing conservative schemes presented in 
this paper we consider as an example the generalised Korteweg-de Vries (gKdV) 
equation 

ut + u xxx + (u^ 1 )^ = 

for an integer p > 3, see for example The case p = 3 is the KdV equation 
(1.1 1, the case p = 4 is known as the modified KdV equation, and p = 6 
is sometimes referred to as the mass critical generalised KdV equation. The 

d 



gKdV can be written as (2.1) with 



H[u] = [ (\ul-\i 
Ju V 2 P 



dx, V = 



dx ' 



The AVF discrete variational derivative (3.6) gives rise to the fully implicit 
scheme (3.5) 



jjn+l _ jjr. 

At 



jjn+1 + y- 




([7«+i)p-i-^({7»)i =0 . 



(6.1) 



After applying the polarising procedure of Section 4.1 to H = J n Gdx we 
get H — J„ G dx which depends on k — \p/2~\ arguments 



G[W!, 



,w k \ 



pk 



I n fc 

. p 



h n- =1 ^ E? =1 i , podd, 



p even. 



After finding the AVF discrete variational derivative from (4.12) we get the 
linearly implicit PAVF scheme (4.13) 

jjn+k _ jjn jjn+k i jjn 



kAt 



(n^c^) 2 ) (if i 1 T ^ rL + i)] x = °> p odd 



nti(^" +j ) 2 (c/™ +fc +c/™) 



(6.2) 



p even. 



Notice that is indeed only appearing as linear terms in this scheme. The 

schemes (1.2) and (1.3) arc found by settings = 3 (fe = 2) in (6.1) and (6.2), 
respectively. Following the procedure of Section [5] one can get a fully discretised 
scheme by replacing U by U and the first and third derivative operators by 5^ 
and 5^ respectively. 

In the Figures 6.1 and 6.2 we compare the conservative methods (1.2) and 



(1.3) with the fully implicit midpoint method 

jjn+1 _ jjn 




= 



and a naive linearly implicit method 



rrn+l _ Tjn jjn+1 , jjn 

_j_ ^ xxx 1 ^ xxx , /jjnjj-n+l\ 

At 2 1 



(6.3) 



(6.4) 
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We test the four methods on a traveling wave solution 



- cb) = ^scch 2 ( ^^{x - ct)\ , c> 



— sech 2 f - — -(t. — 
2 V 2 

using the parameters c = 1, x = (—5,5), Aa; = || and At = 0.1. As an 
indication of the long time behaviour of the presented schemes we consider to 
which extent the methods are able to preserve the shape and propagation speed 
of a traveling wave solution. We define the two quantities 

e shapc = min||(7"-$(--T)|| 2 (6.5) 

T 

and 

£dista„cc = |argmin||[/" - $(• - r)|| 2 - ct n \. (6.6) 

r 

Thus e s hapc measures the shape error of the numerical solution, and Edistance 
measures the error in the travelled distance of the numerical solution. 



We see in Figure 6.1 the fully implicit schemes preserves the shape better 
than the linearly implicit ones, and that the conservative schemes perform better 
than the non-conservative ones. In Figure [672] we see that the linearly implicit 



schemes have a more accurate phase speed than the fully implicit ones. Figure 



6.3 shows the global error as a function of the time step. As expected the plot 
shows that the four methods are second order and that the linearly implicit 
schemes arc inaccurate for large At. In conclusion we see that the linearly 
implicit conservative scheme performs comparably to the other methods while 
being more efficient (the latter is shown in Figure |1.1[ ) . 
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